/*=========================================================================
Elite College as Engines of Upward Mobility: Evidence from Colombia's Ser Pilo Paga
Authors: Juliana Londoño-Vélez, Catherine Rodriguez, Fabio Sánchez
and Luis Esteban Álvarez-Arango

Creation date: June 6th, 2025
--------------------------------------------------------------------------
Figures 7a, 7b, 8a and 8b
=========================================================================*/

use "${data}/data_RD", clear
keep if icfes_per==20142 | icfes_per==20132 | icfes_per==20122 
rename codigo_programa codigo_prog
merge m:1 codigo_prog using "${Intermediate}/va_prog", nogen keep(1 3) keepusing(score_saberpro_5_c_va mw_9_c_va)
replace score_saberpro_5_c_va=. if score_saberpro_5==.
replace mw_9_c_va=. if mw_9_c_va==.

fasterxtile wage0_9_ptile13=wage0_9 if icfes_per==20132, nq(100)
gen wage0_9_top113=(wage0_9_ptile13==100) if icfes_per==20132
fasterxtile wage0_9_ptile14=wage0_9 if icfes_per==20142, nq(100)
gen wage0_9_top114=(wage0_9_ptile14==100) if icfes_per==20142
gen wage0_9_top1=.
replace wage0_9_top1=wage0_9_top113 if icfes_per==20132
replace wage0_9_top1=wage0_9_top114 if icfes_per==20142
drop wage0_9_ptile* wage0_9_top11*

cap drop mean_include*
cap drop cut_saber11 
cap drop cut_sisben
gegen cut_saber11=cut(running_saber11), at(-40(2)40)
gegen mean_include_saber11 = tag(cut_saber11 icfes_per) if eligible_sisben==1
gegen mean_exclude_saber11 = tag(cut_saber11 icfes_per) if eligible_sisben==0

glo cond_saber11	"eligible_sisben==1"
glo cond_sisben		"eligible_saber11==1"
glo no_cond_saber11	"eligible_sisben==0"
glo no_cond_sisben	"eligible_saber11==0"
glo legend_saber11	"SISBEN"
glo legend_sisben	"SABER11"
glo lim_saber11		40
glo lim_sisben		20

gegen cut_sisben=cut(running_sisben), at(-20(1)20)
gegen mean_include_sisben = tag(cut_sisben icfes_per) if eligible_saber11==1
gegen mean_exclude_sisben = tag(cut_sisben icfes_per) if eligible_saber11==0

cap drop bin1 bin2 bin3 bin4
gegen bin1 =  mean(mw_9) if ${cond_saber11} & icfes_per==20142 , by(cut_saber11) // eligible 2014
gegen bin2 =  mean(mw_9) if ${cond_saber11} & icfes_per==20132, by(cut_saber11) // eligible 2013
gegen bin3 =  mean(mw_9) if ${no_cond_saber11} & icfes_per==20142, by(cut_saber11) // non eligible 2014 
gegen bin4 =  mean(mw_9) if ${no_cond_saber11} & icfes_per==20132, by(cut_saber11) // non eligible 2013

#delimit ;
twoway 	(lfit mw_9 running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 < 0 & icfes_per==20132, sort lcolor(gs12))
		(lfit mw_9 running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 >= 0 & icfes_per==20132, sort lcolor(gs12))
		(scatter bin2 cut_saber11 if mean_include_saber11==1 & icfes_per==20132, mc(gs12) ms(Oh) msize(small))
		(lfit mw_9 running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 < 0 & icfes_per==20132, sort lcolor(gs8) lpattern(dash))
		(lfit mw_9 running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 >= 0 & icfes_per==20132, sort lcolor(gs8) lpattern(dash))
		(scatter bin4 cut_saber11 if mean_exclude_saber11==1 & icfes_per==20132, mc(gs8) ms(Dh) msize(small))
		(lfit mw_9 running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 < 0 & icfes_per==20142, sort lcolor(gs4) lpattern(dash))
		(lfit mw_9 running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 >= 0 & icfes_per==20142, sort lcolor(gs4) lpattern(dash))
		(scatter bin3 cut_saber11 if mean_exclude_saber11==1 & icfes_per==20142, mc(gs4) ms(S) msize(small))
		(lfit mw_9 running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 < 0 & icfes_per==20142, sort lcolor(gs0)) 
		(lfit mw_9 running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 >= 0 & icfes_per==20142, sort lcolor(gs0))
		(scatter bin1 cut_saber11 if mean_include_saber11==1 & icfes_per==20142, mc(gs0) ms(T) msize(small)),
		legend(order(3 "Low SES pre" 6 "High SES pre" 12 "Low SES post" 9 "High SES post") position(10) ring(0) col(1)) 
		xline(0, lcolor(red) lpattern(dash)) xsize(8) xlab(-40(20)40) xtitle(Distance to SABER 11 cutoff) ylab(0(0.2)2) graphregion(color(white)) ytitle("")
;
#delimit cr

graph export "${outcomes}/Fig7A.pdf", replace width(20)

cap drop bin1 bin2 bin3 bin4
gegen bin1 =  mean(wage0_9_top1) if ${cond_saber11} & icfes_per==20142 , by(cut_saber11) // eligible 2014
gegen bin2 =  mean(wage0_9_top1) if ${cond_saber11} & icfes_per==20132, by(cut_saber11) // eligible 2013
gegen bin3 =  mean(wage0_9_top1) if ${no_cond_saber11} & icfes_per==20142, by(cut_saber11) // non eligible 2014 
gegen bin4 =  mean(wage0_9_top1) if ${no_cond_saber11} & icfes_per==20132, by(cut_saber11) // non eligible 2013

#delimit ;
twoway 	(lfit wage0_9_top1 running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 < 0 & icfes_per==20132, sort lcolor(gs12))
		(lfit wage0_9_top1 running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 >= 0 & icfes_per==20132, sort lcolor(gs12))
		(scatter bin2 cut_saber11 if mean_include_saber11==1 & icfes_per==20132, mc(gs12) ms(Oh) msize(small))
		(lfit wage0_9_top1 running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 < 0 & icfes_per==20132, sort lcolor(gs8) lpattern(dash))
		(lfit wage0_9_top1 running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 >= 0 & icfes_per==20132, sort lcolor(gs8) lpattern(dash))
		(scatter bin4 cut_saber11 if mean_exclude_saber11==1 & icfes_per==20132, mc(gs8) ms(Dh) msize(small))
		(lfit wage0_9_top1 running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 < 0 & icfes_per==20142, sort lcolor(gs4) lpattern(dash))
		(lfit wage0_9_top1 running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 >= 0 & icfes_per==20142, sort lcolor(gs4) lpattern(dash))
		(scatter bin3 cut_saber11 if mean_exclude_saber11==1 & icfes_per==20142, mc(gs4) ms(S) msize(small))
		(lfit wage0_9_top1 running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 < 0 & icfes_per==20142, sort lcolor(gs0)) 
		(lfit wage0_9_top1 running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 >= 0 & icfes_per==20142, sort lcolor(gs0))
		(scatter bin1 cut_saber11 if mean_include_saber11==1 & icfes_per==20142, mc(gs0) ms(T) msize(small)),
		legend(order(3 "Low SES pre" 6 "High SES pre" 12 "Low SES post" 9 "High SES post") position(10) ring(0) col(1)) 
		xline(0, lcolor(red) lpattern(dash)) xsize(8) xlab(-40(20)40) xtitle(Distance to SABER 11 cutoff) graphregion(color(white)) ytitle("")
;
#delimit cr

graph export "${outcomes}/Fig7B.pdf", replace width(20)

replace icfes_per=20132 if icfes_per==20122
cap drop bin1 bin2 bin3 bin4
gegen bin1 =  mean(score_saberpro_5_c_va) if ${cond_saber11} & icfes_per==20142 , by(cut_saber11) // eligible 2014
gegen bin2 =  mean(score_saberpro_5_c_va) if ${cond_saber11} & icfes_per==20132, by(cut_saber11) // eligible 2013
gegen bin3 =  mean(score_saberpro_5_c_va) if ${no_cond_saber11} & icfes_per==20142, by(cut_saber11) // non eligible 2014 
gegen bin4 =  mean(score_saberpro_5_c_va) if ${no_cond_saber11} & icfes_per==20132, by(cut_saber11) // non eligible 2013

#delimit ;
twoway 	(lfit score_saberpro_5_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 < 0 & icfes_per==20132, sort lcolor(gs12))
		(lfit score_saberpro_5_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 >= 0 & icfes_per==20132, sort lcolor(gs12))
		(scatter bin2 cut_saber11 if mean_include_saber11==1 & icfes_per==20132, mc(gs12) ms(Oh) msize(small))
		(lfit score_saberpro_5_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 < 0 & icfes_per==20132, sort lcolor(gs8) lpattern(dash))
		(lfit score_saberpro_5_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 >= 0 & icfes_per==20132, sort lcolor(gs8) lpattern(dash))
		(scatter bin4 cut_saber11 if mean_exclude_saber11==1 & icfes_per==20132, mc(gs8) ms(Dh) msize(small))
		(lfit score_saberpro_5_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 < 0 & icfes_per==20142, sort lcolor(gs4) lpattern(dash))
		(lfit score_saberpro_5_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 >= 0 & icfes_per==20142, sort lcolor(gs4) lpattern(dash))
		(scatter bin3 cut_saber11 if mean_exclude_saber11==1 & icfes_per==20142, mc(gs4) ms(S) msize(small))
		(lfit score_saberpro_5_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 < 0 & icfes_per==20142, sort lcolor(gs0)) 
		(lfit score_saberpro_5_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 >= 0 & icfes_per==20142, sort lcolor(gs0))
		(scatter bin1 cut_saber11 if mean_include_saber11==1 & icfes_per==20142, mc(gs0) ms(T) msize(small)),
		legend(order(3 "Low SES pre" 6 "High SES pre" 12 "Low SES post" 9 "High SES post") position(10) ring(0) col(1)) 
		xline(0, lcolor(red) lpattern(dash)) xsize(8) xlab(-40(20)40) xtitle(Distance to SABER 11 cutoff) graphregion(color(white)) ytitle("")
;
#delimit cr

graph export "${outcomes}/Fig8A.pdf", replace width(20)

cap drop bin1 bin2 bin3 bin4
gegen bin1 =  mean(mw_9_c_va) if ${cond_saber11} & icfes_per==20142 , by(cut_saber11) // eligible 2014
gegen bin2 =  mean(mw_9_c_va) if ${cond_saber11} & icfes_per==20132, by(cut_saber11) // eligible 2013
gegen bin3 =  mean(mw_9_c_va) if ${no_cond_saber11} & icfes_per==20142, by(cut_saber11) // non eligible 2014 
gegen bin4 =  mean(mw_9_c_va) if ${no_cond_saber11} & icfes_per==20132, by(cut_saber11) // non eligible 2013

#delimit ;
twoway 	(lfit mw_9_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 < 0 & icfes_per==20132, sort lcolor(gs12))
		(lfit mw_9_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 >= 0 & icfes_per==20132, sort lcolor(gs12))
		(scatter bin2 cut_saber11 if mean_include_saber11==1 & icfes_per==20132, mc(gs12) ms(Oh) msize(small))
		(lfit mw_9_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 < 0 & icfes_per==20132, sort lcolor(gs8) lpattern(dash))
		(lfit mw_9_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 >= 0 & icfes_per==20132, sort lcolor(gs8) lpattern(dash))
		(scatter bin4 cut_saber11 if mean_exclude_saber11==1 & icfes_per==20132, mc(gs8) ms(Dh) msize(small))
		(lfit mw_9_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 < 0 & icfes_per==20142, sort lcolor(gs4) lpattern(dash))
		(lfit mw_9_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${no_cond_saber11} & running_saber11 >= 0 & icfes_per==20142, sort lcolor(gs4) lpattern(dash))
		(scatter bin3 cut_saber11 if mean_exclude_saber11==1 & icfes_per==20142, mc(gs4) ms(S) msize(small))
		(lfit mw_9_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 < 0 & icfes_per==20142, sort lcolor(gs0)) 
		(lfit mw_9_c_va running_saber11 if abs(running_saber11)<${lim_saber11} & ${cond_saber11} & running_saber11 >= 0 & icfes_per==20142, sort lcolor(gs0))
		(scatter bin1 cut_saber11 if mean_include_saber11==1 & icfes_per==20142, mc(gs0) ms(T) msize(small)),
		legend(order(3 "Low SES pre" 6 "High SES pre" 12 "Low SES post" 9 "High SES post") position(10) ring(0) col(1)) 
		xline(0, lcolor(red) lpattern(dash)) xsize(8) xlab(-40(20)40) xtitle(Distance to SABER 11 cutoff) graphregion(color(white)) ytitle("")
;
#delimit cr

graph export "${outcomes}/Fig8B.pdf", replace width(20)

